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Abstract 

In this paper we discuss some features of the BCRE model. We show that this 
model can be understood as a mapping from a two-dimensional to a one-dimensional 
problem, if some conditions are satisfied. We propose some modifications that (a) 
guarantee mass conservation in the model (what is not assured in its original form) 
and (b) correct undesired behaviors that appear when there are irregularities in the 
surface of the static phase. We also show that a similar model can be deduced both 
from the principle of mass conservation (first equation) and a simple thermodynamic 
model (from which the exchange equation can be obtained). Finally, we solve the 
model numerically, using different velocity profiles and studying the influence of the 
different parameters present in this model. 
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1 Introduction 



Since the end of 1980s, when Bak et al. used a sandpile as a paradigm for 
self-organized criticality[6], the interest in granular materials experimented 
a revival. The topic, however, is not recent. The first studies with granular 
materials date back to 1773, when Coulomb first observed that this kind of 
matter could stand in equilibrium in piles at certain specific angles and, after 
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that, other famous physicists studied the topic. Faraday discovered the con- 
vective instability in vibrating grains and Reynolds introduced the notion of 
dilatancy, just to cite some examples. 

The study of granular materials is important since their applications in in- 
dustry is wide enough to cover areas as distinct as civil construction and food 
transportation. The study of flows can uncover the behavior of dunes, sand- 
storms and avalanches. But those are just a few examples. Granular materials 
are present everywhere in nature, and in several branches of industry. Chem- 
ical industries, pharmaceutics, mining, geology are just some more examples 
of other areas where the study of grains can lead to important results. 

In 1994, a paper by Bouchaud, Cates, Ravi and Edwards [1] presented an 
important model that became known as the BCRE model. This model was 
successful in describing the qualitative behavior of flowing grains, and had the 
additional advantage of being very simple. It assumes that a two-dimensional 
sandpile (figure 1), with rolling grains on its surface, can be divided in two 
"phases" (a static phase h and a rolling phase R) and propose two coupled 
partial differential equations to model their behavior: 



dR(x,t) 
di 



D 



dR(x,t) 
dx 



+ r(R,k), 



(i) 



r (r, h) = -r (x, t) 



dh(x,t) d 2 h(x,t) 

7 TT^ + K 



dx 



dx 2 



dh (x, t) 



(2) 



where 7, n > and h = h + :rtan# r , 6 r is the angle of repose. 

The first equation defines how the profile of the rolling phase evolves in time, 
and the second equation determines the profile of the static phase by setting 
the form of the exchange between rolling and static grains, depending on the 
local slope of the pile. The model is phenomenological, and is obtained from 
educated guesses about the properties observed in real sandpiles and desired 
for real grains. 

This model gave important clues to how we could describe some interesting 
phenomena occurring in granular flow. A variety of papers utilized these equa- 
tions to model the behavior of avalanches, stratification and flows in general [3- 
5,9,10]. 

The original model, however, is very simplified and need some improvements 
to describe more general situations. Some problems of consistency must also be 
solved if we want to have a model that do not violate some general principles, 
as mass conservation. The model, also, still lacks from a derivation from first 
principles or from a microscopic point of view. 
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In this paper we will address some of the points mentioned above. In section 2, 
we will show that an equation equivalent to equation (1) can be obtained from 
the principle of mass conservation, under the assumption that the densities 
of the static and rolling phases are constant in the vertical coordinate. In 
section 3, we discuss some limitations of equation (2) and present some possible 
alternatives to it. In section 4, we analyze the consequences of considering 
different velocity profiles for the rolling phase and discuss the relation of them 
with other models for granular materials present in literature, evaluating the 
role played by some of the parameters present in the model. In section 5, we 
present a new and simple model to describe the mechanism that underline the 
exchange of grains between static and rolling phases. From it we were able to 
deduce an equation of the type of equation (2). Finally, in the last section, we 
summarize our conclusions. 



2 Mass Conservation 

Consider, for instance, a two-dimensional sandpile with a rolling and a static 
phase. The rolling phase is located above the static phase, and slides over it 
(see figure 1). Let us assume that the sandpile can be treated as a continuous 
medium, p r being the density in area of the rolling phase. The mass inside an 
interval x Q to x Q + 5x, where 5x is small, is then given by 




with 



h (x, t) = height of the static phase in a point x at time t, 
r (x, t) = height of the rolling phase in a point a; at a time t. 



We define the linear density of the rolling phase by 




Remembering that, in general, 





X 
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we can write (3) as 



h(x ,t)+r(x ,t) 

R(x ,t)= J p r (x ,y,t)dy. (5) 

h{x ,i) 



When p r is independent of the vertical coordinate y, equation (5) becomes 

R(x,t) = p r (x,t) r (x,t) . (6) 

Repeating this procedure for the static phase, we obtain 

h(x a ,t) 

S (x , t) = hm = J p s (x , y, t) dy, (7) 



o 



and 



S(x,t) =p s (x,t) h(x,t), (8) 



where m s is the mass and p s (again independent of y) is the density in area 
of the static phase. 

Equations (6) and (8) define a one-to-one mapping from h and r to S and 
R respectively, mapping the two-dimensional sandpile in a one-dimensional 
problem in the case that the densities of the two phases do not vary with y. 
A particular case that follows in this category, used frequently in literature, 
is when p r — p s — constant. The possibility of this mapping is not really 
unexpected, since, with this properties, the variables involved will only depend 
on the horizontal coordinate x and time. 

Under the assumption of a continuous model, we can write, for the rolling 
phase, a mass conservation equation. For a one-dimensional fluid of density R 
that flows under a velocity field v (if v also does not vary in the y coordinate) 

dR 9 , m ^ 

m + d-x {vR) = Q > 



where Q represents the sources or sinks of this fluid. The BCRE model allows 
the exchange between rolling and static grains. So, when describing the rolling 
phase, the static phase acts like a source/sink at every point (the extra mass 
gained by the rolling phase equals the mass lost by the static phase). We can 
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then write 



(9) 



Comparing (1) and (9), we can immediately see that the last is a slightly 
modified version of the first, where, instead of the height h, we are working 
with the density S and the diffusion term (if there is one), will now depend 
on the specific shape of the velocity field v. 

Equation (9) is better than equation (1), since it guarantees mass conservation 
for any form of v (what does not happen with equation (1)). If we assume, 
for instance, a constant velocity field in (1) (as has already been done in some 
previous works [[9,3,4]), the diffusing term must be dropped out in order to 
ensure mass conservation. 

In 1997 Makse [9] applied the BCRE model to a mixture of two grains. He 
wrote equation (1) — with constant velocity and without diffusion — for each 
kind of grain 

dRj dRj _ _ „ , _ „ . 

_ + t ,,_ = r„ . = 1,2, (10) 



and simulated this model with v\ — v 2 , showing that, depending on the value 
of some parameters, the grains would segregate or stratify. The assumption 
V\ = v 2 was not justified in that paper. To exemplify the advantages of our 
equations, applying (9) to a mixture of two grains, and considering that the 
mixture is rolling with a constant average velocity v, we obtain 

d(R 1 + R 2 ) , d 

— m — + «^(*i + ^) = r, 



that is similar to the equations obtained by Markse (10). But now it is easy 
to see that v\ = v 2 is not an ad hoc assumption, but rather, an imposition of 
the model. 

One more advantage of (9) is that now we can obtain a whole series of different 
granular flow regimes by varying the velocity profile of the rolling phase. We 
will discuss this point better in section 4. 



3 Exchange of grains between static and rolling phases 



We will now make some considerations about the second equation of the BCRE 
model, that we call exchange equation. In analogy with (9) we will first write 
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this equation (2) in terms of the new variables R and S 



(dS „\ d 2 S 

7 l — + Pa tme r 



In [5], Boutreux proposed a modified version of this equation, to take into 
account a shielding effect, present on the upper grains of the rolling phase, 
due to the lower grains of this same layer. Adapted to the variables R and S, 
this shielding effect can be introduced in equation (11) as 



dS_ 
~dt 



R+i> 



7 I + pst&nOr 



+ K 



d 2 S 
dx 2 



;i2) 



where £' is a small constant related to the thickness of the layer of rolling 
grains that indeed interact with the static phase. Note that, if it! ~ £' (a thin 
rolling phase), R£'/(R + f ) ~ R, but if R > then B£'/(R + CO ~ 

Equation (12) is still not adequate to describe what happens close to the 
interface, if there are irregularities on it. As the grains flow, they can erode 
part of the static phase, creating a (sometimes big) crater with a positive slope 
in the right border (see figure 2). To see that, let us analyze a particular case, 
where the densities are constant. 

The term dS/dx is directly related to the slope of the pile (given by dh/dx). 
If it is negative, the pile is inclined to the right (and if it is positive the pile 
is inclined to the left). There is no problem when the slope is negative. As 
expected, for a local slope above the repose angle, we have erosion (and for a 
slope below it we have acrescion). But, when this term is positive, we always 
have acrescion, what is not a reasonable behavior. To correct that we suggest 
a modification, where we consider the minus sign of the modulus of dS/ dx 



dS_ 
~dt 



R + ? 



7 p s tan# r 



8S_ 

dx 



+ K 



d 2 S ~ 

dx 2 



(13) 



4 Diffusion and velocity profiles 



One interesting, important and observed property of the BCRE model is that 
it includes a diffusion of the rolling grains. The presence of diffusion in this 
phase is observed and, in fact, very obvious. This diffusion, however, leads to 
some constraints on the velocity profile v. 

The correct expression for the velocity field v, in equation (9), should, in prin- 
ciple, be derived from a momentum conservation equation. However, writing 
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such an equation is not an easy task, since it should take into account all the 
interactions between the grains and would depend on the stress tensor of the 
material. In general, we can only say that the velocity field is a function of 
x and t that depends on a variety of factors. Most of the works present in 
literature have assumed a constant velocity profile, for simplicity. 

We will, in this section, analyze two possible functional forms for v , showing 
some results of computer simulations and studying its effects in the shape of 
the pile. 

In all computer simulations we integrate the BCRE equations numerically by 
means of a finite difference scheme (see appendix) and constructed an online 
animation of the profile of both phases in real time. The figures presented bel- 
low are screenshots of the animation generated by the program. We considered 
the following profiles: 

(1) v = v(R) = ad x R + PR + 5 , where a, j3, 5 are constants. 

Here v is considered as being a function of R only, and a kind of gradient 
expansion was done. Note that, in this case, v is an implicit function of x and 
t (v = v(R) where R = R(x,t)). Substituting in (9) we obtain 

OR d . _ dR (dR\ 2 , _„ dR B 2 R 
— + — (vR) = — + a — + (2(3R + 5) — + aR—r. (14) 



dt dx dt \dx ) dx dx 



Note that this functional form of v includes a diffusion term. It can also be 
seen from (14) that: 

i) The diffusion coefficient is proportional to R. This means that the higher 
the pile, the more it will diffuse. This is reasonable, and can be understood as 
a consequence of gravity. 

ii) Equation (14) presents a non-linear term with a constant coefficient a, 
(a ^ 0). When diffusion is small and can be neglected (as in most cases studied 
in literature), the non-linear term is also unimportant, and the equation above 
is in agreement with previous works. Diffusion, however, affects the profile of 
the pile. Figure 3 shows the effects of this term. It compares the profile of 
a pile of grains, at the same time and for the same initial condition, for: 
(a) the original model with constant velocity and (b), (c) equation (14) with 
v = ad x R, for two different values of a. We can see that the final shape of the 
bump is quite different if a is not too small. 

iii) This equation has also an advective term, i.e., a term on the first derivative 
of R with respect to x. Its coefficient depends linearly on R (for (3 ^ 0) 

3 An advective term like that had already been proposed in other papers, as in [7] 
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Figure 4 shows how it affects the profile of the pile. We can see that there is 
a tendency to the formation of shock fronts in one of the sides of the bump. 



(2)v = f(d x h) 

In a real pile of grains, there may be irregularities with positive slope, due to 
erosion. In this case, the velocity must depend on the sign of the slope, or else 
we will have avalanches climbing up the pile at the points with positive slope, 
with the same velocity as in the negative slope side. To correct this defect, we 
considered the following form for v 



v = < 



ad x R + pR + 5, ifd x h<0 
v s , if d x h > 



(15) 



where v s is a constant. 

For d x h < this expression is equivalent to the previous form of v. But, for 
d x h > 0, the rolling grains, now, meet a barrier of static grains, and move up 
this barrier with a constant velocity. There is no diffusion in this case, since 
the velocity is constant. But now this is a desired property: the grains slowly 
accumulate in the barrier (and do not diffuse). Figure 5 shows the profile of 
the pile at the same time when (a) the velocity is independent of the slope 
and (b), (c) v depends on the slope according to (15) for two different values 
of v s . The static phase is gray, and has been settled to an irregular shape to 
amplify the consequences of the modifications introduced. 

When we were at the end of this work, we came to know about a work of 
Herrmann and Sauermann, on the behavior of the Barchan Dunes [2], that 
was also based on the BCRE model, and dealt with some of the points we 
present in this paper. They present a two-dimensional version of (9) but they 
do not deduce it. In particular, there is no mention of the need of independence 
of the densities with the y coordinate. Indeed, we believe that, (although it is 
not explicitly written) the simulations were performed under constant densities 
which is a particular case of y independence. They also used a slightly different 
version of (2), where the modulus of the slope of the static phase were also 
employed. 



5 Simple Model for the exchange of grains between rolling and 
static phases 



It is possible, on the bases of a naive model, to deduce an equation to describe 
the exchange of grains between rolling and static phases. Remember that 
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the densities of the rolling and static phases are different, the latter being 
higher than the former. As the rolling phase rolls over the static phase, the 
friction between them takes energy out of the rolling grains. Part of this energy 
becomes heat, but the rest of it is transferred to the grains of the static phase. 
This energy will agitate these grains and the density of the static phase right 
below the interface will increase until it attains the critical dilatancy and starts 
to move, becoming part of the rolling phase. This process is similar to a solid 
to liquid first order phase transition, the static phase playing the role of the 
solid (receives energy and starts to "melt"). Indeed, there are experimental 
evidences that, at least in the case when the transition is induced by tilting, 
it does have features of a first order phase transition [8]. We will assume that 
this analogy is valid. We can than calculate the mass of static grains that will 
"melt" (that is, receives energy and starts to roll) using an analogy with the 
latent heat equation 

SQ = L Sm s , 

where SQ is the energy acquired from the rolling phase, Sm s is the mass of 
static grains that melts and L is a constant, analogue to the latent heat. 

Let us now focus on what happens in a little interval Sx of the horizontal 
coordinate of the pile (see figure 1). The mass that is melted is a portion 
of the total mass of static grains. However, not all the static grains receive 
energy from the rolling phase. The upper grains shield the lower grains from 
the contact with the rolling phase. We will consider that the mass that can 
actually receive energy (and, therefore, melt) is a fraction £h of the static 
phase (we will justify this assumption, now adopted for the sake of simplicity, 
ahead) . 

So, the mass that can be melted is given by 

S 

m s = p s V = £hSx p s — £ h Sx — , 

h 

where equation (8) were used. Then, Sm s = SS £Sx and 

5Q = L£5SSx, (16) 

where SS is the variation in S caused by the melting of the static phase. Note 
that when the pile is too high, internal forces and gravity act to de-estabilize 
it, increasing the melting. This can justify the assumption, made above, that 
m s is proportional to h (instead of being a constant layer). 

If all the energy to melt the static phase comes from the rolling phase, and 
if it is a fraction of the kinetic energy lost due to friction in the interface, we 
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have 

5Q = c5K, (17) 

where SK is the kinetic energy lost by the rolling phase. But K = p 2 /2m r , 
where m r is the mass of the rolling phase and p its momentum. So we have 

5K= 2pSp = (m r v)6p =v 
2m r m r 

Remembering that m r = RSx, supposing that the rolling phase transfers 
energy to the static phase only by friction, and that the friction is proportional 
to the weight of the rolling phase at x, we have 

dp 

— — /im r g =>- 5p — fj,m r g St — /i g R5x St, (19) 
at 

and, from (16), (17), (18) and (19), we have 

L £ SS Sx = c n gv RSt Sx. 

That, in the limit St, SS — >• 0, gives 

-^ = (3vR, (20) 

where p = ——— is a constant. 

This is a quite simple expression for the exchange equation. It can assume a 
variety of forms, depending on the velocity field v. It incorporates the velocity 
field explicitly, indicating that the exchange of grains between both phases 
depends on the exact shape of v, what is very reasonable, since the velocity of 
the rolling grains interfere directly with the energy lost in the collisions, that, 
ultimately, are responsible for the transformation of static grains in rolling 
grains. 

Note also that if (20) is substituted in the mass conservation equation, we get 

dR dR _ 

dt V dx ^' 



where 



>-=)* 
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When the velocity is a function of x, t and R only, q will be, in principle, also 
a function of these three variables (because R is a function of x and t only), 
and the resulting equation will be a well known quasi-linear partial differential 
equation of first order in R, that can be solved by the method of characteristics 
given by the simple system 

dt dx dR 

1 v(x,t,R) q(x,t,R) 



The only difficulty is that, if v has an explicit dependence on R, the variable q 
will have a dependence in dR/ dx, what may make the system of characteristics 
difficult to be solved analytically in practice. However, a solution for R(x, t) 
will give consequently a solution for S(x,t) by means of (20). We intend to 
further explore this point in a following paper. 



Furthermore, if equations (20) and (13) are equivalent, the velocity field must 
assume the form of 



v = 



e/P 

R + ? 



7 p s tan 9 r 



dS_ 

dx 



+ K 



dx 2 



This velocity field has some interesting features. First, note that v depends on 
the slope of the pile. The first term inside the square brackets becomes positive 
for a slope below the angle of repose and negative above it, what means that 
the velocity is higher where the pile is steeper. Second, the factor (R + £')~ 1 
suggests that v is inversely proportional to the weight of the pile of rolling 
grains, what is reasonable. 



6 Conclusions 



In the first two sections of this paper, we pointed out some problems with the 
equations of the BCRE model, suggesting the following modification: 



dR d , nx dS 

-dl + Yx {vR) = ~W 



(21) 



dS_ 
~dt 



R + ? 



7 p s tan# r 



dS_ 

dx 



K 



d 2 S ~ 

dx 2 



(22) 



The first equation explicitly assures mass conservation. The second equation 
was modified to take into account the signal of the slope in the static phase. 
In addition, we introduced the new variables R and S, giving a precision 
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definition for them, a point that was not totally clear in the literature until 
now. 

We used the corrected equations (21) and (22) to simulate two possible velocity 
fields, and found acceptable results. 

Finally, we have proposed a model for the exchange of grains between the 

dS 

rolling and static phases. From it, we obtained the alternative — — = 3v R for 

at 

the exchange equation, that is simpler, includes explicitly the velocity field, 
and lead to interesting results. 

The model presented is very simple, and we are aware that many of the hy- 
pothesis made must be examined in more detail. We assumed that the friction 
is proportional to the weight of the rolling phase, and this is surely oversim- 
plified. Probably there is a more complicated dependence on other parameters 
of the model as well. For instance, it is reasonable to suppose that the energy 
transfered to the static phase also depends on the shape of the grains, its 
density, the toughness of the material and a variety of other factors. 

Also, we think that the fraction of static grains that receives energy from 
the rolling phase is probably a more general function of h, not to speak of a 
dependence on the other variables of the model. 

We have also neglected the inverse process, i.e., the transformation of rolling 
grains in static ones. Our model may be good to describe an avalanching 
process, where the inertia of the rolling grains is large, but may fail to describe 
a more general situation. 

We hope to address this points in a following paper. However, we think it is 
amazing and interesting that a somewhat rich expression for v, like (20), can 
be obtained from such a naive model. 

Acknowledgements: R. C. A. acknowledges the financial support of the Brazil- 
ian agency FAPESP. 



7 Appendix 

The equations of the BCRE model were integrated with the operator splitting 
method [11]. Suppose a differential equation of the form 

du 
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where L(u) = J2l=i Li(u) is a generic non-linear operator that can be written 
as a sum of m other operators, and u is a function of x and t. If we have a 

good method to integrate each of the equations — = Li(u), then u n+1 can be 
obtained through m successive time steps: 



u n+(l/m) = ^ 
M n+(2/ m) = L2 ^n+(l/ m)) ft) 



u n+1 



L m {y^Hm-l)/m^ ft 



For the first equation of the model (and its variants), the Li operators are of 
the form 



dr <9 2 r / dr \ 2 

L 1 (r)=f(r) — , L 2 {r)=g(r)—, L 3 (r) = q(r) and L 4 (r) = A;l — I, 



where /, (7 and g are arbitrary functions of r , and k is a constant. 

L\ and L2 where integrated with variants of the Crank-Nicholson method; the 
operator L 3 was integrated with a fourth order Runge-Kutta procedure and 
finally L 4 was integrated with a FTCS finite-difference scheme. 

The second equation can be split into two operators of the form: 



U (h) = a|£, L 2 (h) = bf x , 

where a and b are constants. Now L\ was integrated with the aid of Crank- 
Nicholson and the operator L 2 with the aid of a two-step Lax-Wendroff [11] 
procedure. 

The program was run on line in a Intel Pentium III 400 Hz or in a AMD K6II 
400 Hz (about 10 seconds for a single run in both cases). 
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Figure Captions 



Figure 1 In the BCRE model, a two-dimensional pile with flowing grains 
is divided into two "phases" (static and rolling phases). They are described, 
respectively, by the variables h and R = pr, where h and r are the heights of 
the static and rolling phases, respectively. 

Figure 2 Eventually the flow of rolling grains can erode the static phase 
creating a crater, and generating regions of positive slope. 

Figure 3 Profile, at equivalent time t and same initial conditions, of the 
evolution of a pile in cases where: (a) the original model was considered with 
constant velocity (v = 1); (b) our suggestion mass conservation (equation (9)) 
with a = — 0.01 and (c) a = — 0.04. Gray area corresponds to the static phase 
and white area to the rolling phase. 

Figure 4 Effects of diffusion. Profile, at equivalent time t and same initial 
conditions, of the evolution of a pile with a velocity profile given by v = 
ad x R + (3R + 5, with 5 = 1 , a = 0.01 and (a) p = 2 and (b) (3 = - 1 (A 
negative j3 means that the thicker the rolling phase,the more difficult it is to 
th grains to roll). Gray area corresponds to the static phase and white area 
to the rolling phase. 

Figure 5 Profile, at equivalent time t and same initial conditions, of the 
evolution of a pile in the cases where (a) the velocity is indepedent of the sign 
of slope and (b) the velocity depends on the sign of the slope and v s = 0.9 
and (c) v s = 0.4. Gray area corresponds to the static phase and white area to 
the rolling phase. 
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